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Abstract 

We study a rate-model neural network composed of excitatory and inhibitory neurons 
in which neuronal input-output functions are power laws with a power greater than 1, 
as observed in primary visual cortex. This supralinear input-output function leads to 
supralinear summation of network responses to multiple inputs for weak inputs. We 
show that for stronger inputs, which would drive the excitatory subnetwork to insta- 
bility, the network will dynamically stabilize provided feedback inhibition is sufficiently 
strong. This dynamic stabilization yields a transition from supralinear to sublinear 
summation of network responses to multiple inputs. We compare this to the dynamic 
stabilization in the "balanced network", which yields only linear behavior. We more 
exhaustively analyze the 2-dimensional case of 1 excitatory and 1 inhibitory population. 
We show that in this case dynamic stabilization will occur whenever the determinant of 
the weight matrix is positive and the inhibitory time constant is sufficiently small, and 
analyze the conditions for "supersaturation" , or decrease of firing rates with increasing 
stimulus contrast (which represents increasing input firing rates). In work to be pre- 
sented elsewhere, we show that this transition from supralinear to sublinear summation 
can explain a wide variety of nonlinearities in cerebral cortical processing. 
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1 Introduction 



We have recently found, in work that is yet unpubhshed except as abstracts (Miher and Rubin 2010, 
2011, Rubin and Miller 2010, 2011) that a large set of response properties of cells in primary visual 
cortex (VI) and other sensory cortical areas can be understood from a very simple circuit motif. 
The response properties have in common a change in integration with increasing input strength, so 
that responses to weak inputs sum supralinearly while those to stronger inputs sum sublinearly. 

One set of properties involve contextual modulation or "surround suppression" . A visual sensory 
neuron has a classical receptive field (CRF), corresponding to the region in which appropriate visual 
stimuli will drive the neuron's responses. The size of the CRF does not change with input strength 
(Song and Li 2008). Stimuli outside the CRF can modulate responses to CRF stimuli, although 
they cannot drive responses, and typically are suppressive. However, the nature of the surround 
influence can vary with input strength (Polat et al. 1998, Sengpiel et al. 1997). A size tuning curve is 
obtained by centering an effective stimulus on the CRF center and studying response vs. stimulus 
radius. The summation field size is the stimulus size evoking peak response. This summation 
field size shrinks with input strength, as represented by stimulus contrast (Anderson et al. 2001, 
Cavanaugh et al. 2002, Sceniak et al. 1999, Shushruth et al. 2009, Tsui and Pack 2011). This 
means that regions of the surround are changing from facilitating to suppressing with increasing 
input strength. 

Another set of properties involve sublinear summation of the responses to multiple stimuli: the 
response to two simultaneously presented stimuli can be closer to the average than the sum of 
the responses to the stimuli presented individually. We refer to this property as "normalization", 
because it is the most prominent of a set of nonlinear response properties that have been given 
that name (reviewed in Carandini and Heeger 2011). In at least some cases, this summation 
becomes supralinear when inputs are weak (Heuer and Britten 2002, Ohshiro et al. 2011). If one 
thinks of surround suppression as representing the response to simultaneous presentation of a center 
stimulus that normally by itself evokes a certain response and a surround stimulus that normally by 
itself evokes zero response, then surround suppression can be thought of as an example of sublinear 
summation. Similarly, facilitation by the near surround for weak inputs then represents supralinear 
summation. 

We have found that these and other response properties can be understood in some detail from 
a simple model. We consider a network of excitatory (E) and inhibitory (I) neurons, extended 
across a 1-D or 2-D space. The strengths of each type of connection - E^E, E— >I, I^E, I^I - 
fall off as functions of distance. As in biology, I projections are shorter-range than E projections. 
We are guided by previous results showing that the inhibition received by cells is decreased when 
they are being suppressed by a surround stimulus, relative to their response to a CRF stimulus 
alone (Ozeki et al. 2009), and correspondingly showing that the firing of inhibitory cells, like that 
of excitatory cells, is suppressed by surround stimuli (Song and Li 2008). These results led to the 
conclusion that the E— t-E connections must be sufficiently strong that, when the network is being 
driven by a CRF stimulus, they would render the network unstable in the absence of feedback 
inhibition (Ozeki et al. 2009), a conclusion also supported by other work (London et al. 2010). We 
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term such a network an inhibition-stabilized network or ISN. 

We then add to this the fact that individual neurons have a supralinear, power-law input-output 
function. This is based on intracellular recordings in anesthetized cat primary visual cortex (VI) 
showing that a neuron's instantaneous firing rate is well described as a power law function of its 
instantaneous mean voltage relative to rest (rates and voltages measured in 30 ms bins) with powers 
ranging from 2 to 5, and that this holds true over the entire dynamic range of neuronal response 
to visual stimuli (Finn et al. 2007, Priebe and Ferster 2005, 2006, Priebe et al. 2004).^ This power 
law relationship is predicted on theoretical grounds when mean input is subthreshold and spiking is 
driven by input fluctuations (Hansel and van Vreeswijk 2002, Miller and Troyer 2002), as appears 
to be the case in VI (Anderson et al. 2000b). 

We find that this supralinear network can explain a wide variety of response properties including 
those described above. Here we mathematically analyze the model. We focus particularly on 
exposing the origins of the transition in model behavior from supralinear to sublinear summation 
with increasing input strength, which occurs as the excitatory subnetwork becomes unstable and 
is stabilized by feedback inhibition. Hence we refer to the network as the stabilized supralinear 
network or SSN. We also conduct a more detailed analysis of the 2-dimensional case consisting of 
a single excitatory and a single inhibitory population. 



2 Setup: Equations for the Supralinear Network 

We take r = ^ ^ to be the AT-dimensional vector of neuronal firing rates, ordered so that the 

top Ne neurons, represented by rg, are all excitatory neurons, and the remaining Nj neurons, 
represented by r/ are inhibitory neurons, Ne + Nj = N. (We refer to the units in our model as 
"neurons", but, as discussed below, the equations represent average firing rates and so excitatory or 
inhibitory units may be better understood as local interconnected groups of excitatory or inhibitory 
neurons, over which the average is taken.) The matrix of connections between the neurons is 
/ Wee -Wp;/ \ 

W = where Wxy is the matrix of connections from neurons of type Y (E or 

I) to neurons of type X and has non-negative entries. The feedforward input to the neurons in the 
hE 
hi 

We study the simplest standard firing-rate-model equations (reviewed in Ermentrout and Ter- 
man 2010, Chapter 11; Gerstner and Kistler 2002, Chapter 6; Dayan and Abbott 2001, Chapter 



network is h 



'^We are assuming that mean voltage is linear in the input. Nonlinearities such as spike-rate-adaptation currents 
could complicate this picture. Wc also arc ignoring the fact that the power increases with contrast, because the 
noise level decreases with contrast (Finn et al. 2007), which yields increasing powers (Hansel and van Vreeswijk 2002, 
Miller and Troyer 2002). However the picture we describe in this paper primarily concerns stabilization against the 
otherwise explosive nonlinearity of a supralinear input-output function. Thus, the picture should hold so long as the 
input-output function is supralinear over the cell's dynamic range, as expected for fluctuation-driven spiking - the 
closer the cell is to threshold, the greater the increase in spiking driven by a given increment of input. 
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7), in which a neuron's firing rate approaches a nonhnear function of its input with first-order 
dynamics: 

tT— = -r + f(Wr + h) (1) 

Here T is a diagonal matrix of relative time constants, i.e. the time constant of the i*'' neuron is tTh. 
f is a vector function of a vector argument that acts elementwise on its argument, (f (v)), = fi{vi), 
for some scalar functions of a scalar variable, fi, where Vi is the i^^ element of v. These rate 
model equations do not capture fast time scales that arise in spiking networks, and cannot capture 
synchronization of spikes across neurons, but tend to be reliable in describing steady states or 
slower aspects of dynamics when neurons spike asynchronously. We will focus on the steady state 
and its stability. 

We will study the case in which the fi are identical for all elements, fi = f, and / is a rectified 
power law with power n > 1: 

fix) = k{[xur (2) 

where [x]-^- = x, x > 0; = 0, otherwise. We will summarize this by saying 

rT^ = -r + A;(Wr + h)-" (3) 

where v" is the vector with i*'* element ([uj].).)" (the period in the exponent .n, based on Matlab 
notation, is to indicate that the operation is done element-by-element rather than to the vector as 
a whole) . A power- law relation between the mean input and mean response arises in the case that 
spiking is driven by input fluctuations (Hansel and van Vrceswijk 2002, Miller and Troyer 2002), 
and similarly it is observed in VI as the relation between the trial-averaged mean voltage and mean 
response. 

We now change variables to dimensionless ones. This allows us to determine the dimensionless 
combinations of parameters on which model behavior depends and in which expansions for small 
or large values may be undertaken. We let ip = ||W|| where ||W|| is some matrix norm or other 
measure of the size of W, and write W = V'J with J dimensionless and ||J|| = 1. Similarly we 
let c = |h| and write h = eg with g dimensionless and |g| = 1 (again, |g| indicates some measure 
of the size of a vector, e.g. a vector norm). Note that c and tpr have the same units, so that 
i/jr/c is dimensionless, and that fee" and r have the same units, so that kc^/{c/il)) = kc^'^i/j is 
dimensionless. We thus define the dimensionless variable and parameter: 

y = rV'/c (4) 
a = kd^-^i^ (5) 



Then equation 3 becomes 
Thus, given J, g, T, and n, the dynamics depends only on the single parameter a. 



rT^ = -y + a(Jy + g)- (6) 
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The fact that Eq. 6 has a single a for all neurons is quite general: if neuron i had parameter CKj, 
this could be replaced with a by multiplying all weights Jij and inputs gi to neuron i by (^)^''", 
leaving the form of the equation unchanged. However the fact that the equation has a single n for 
all neurons is a real restriction. Consideration of ra's that vary between neurons or between neuron 
types remains a question for future study. 

Note that we can rewrite the input to r in our model as (fcnWr + /cncg) ", so that the effec- 
tive recurrent weights are kn'W and the effective input strength is knc. Then a = kd^^^ip = 

(A;^||W||) (kUy~ , that is, a = (recurrent weight) (feedforward weight)" ^. Note also that 
whether the input is dominated by feedforward input eg or recurrent input Wr is determined 
by the size and structure of y for a given a (because Wr + h = c( Jy + g) , so that the balance 
depends only on the relative sizes of Jy vs. g), and is not impacted at all by the ratio c/ip, which 
naively might be thought to determine the feedforward/recurrent balance. For a given a, this ratio 
simply scales r (r = {c/ilj)y). 

We will focus on the equation for the steady-state: 

y = a(Jy + g)- (7) 

However, in considering stability of the steady state we will need to use the dynamical equation 6. 

3 Scaling Argument 

In this section, we show that the supralinear network generically makes a transition from supralinear 
summation of responses to multiple sets of feedforward inputs for weak inputs (q <C 1) to sublinear 
summation for stronger inputs {a S> 1), with the transition occuring for a of order of magnitude 
1, for which we use the standard notation a ~ 0(1)- This transition occurs because of dynamic 
stabilization by feedback inhibition of an otherwise explosive network. We then compare this 
stabilization to that in the balanced network model of Van Vreeswijk and Sompolinsky (1998). 

3.1 Scaling for small a 

For a ^ 1 we expect the steady state to satisfy y ag ", since then the Jy term is small relative 
to the g term and so adds only a small correction to this solution. More generally, we can write a 
formal expression for the steady state: 

y = a(g + aJ(g + aJ(g + aJ(. . .)-)-)-)- (8) 

or 

r = k{cg + kW{cg + A;W(cg + kW{. . .)■")■")■")■" (9) 

where the ellipses indicate infinite repetition of the pattern. Assuming quantities in the parentheses 
are positive so that we can ignore rectification, which they will be for sufficiently small a, equation 8 
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can be converted into an infinite series in increasing integer powers of a with dominating (lowest- 
order) term ag'".^ The terms multiplying will involve factors of g interspersed with J's, with 
the sum of the powers on the g's equal to p{n — 1) + 1. Similarly for r, one obtains a series involving 
an infinite set of powers of c, with lowest-order term k{cg)''^ and higher-order terms proportional 
to k^c^^'^^^^'^^ and involving a set of g's with summed power also equal to p{n — 1) -|- 1. If this 
series converges, which it will for sufficiently small a, it will give a steady state solution. 

Thus, for small a, feedforward inputs sum supralinearly to produce responses. Intuitively: the 
effective connection between two neurons tells how much the steady-state postsynaptic rate changes 
for a given change in steady-state presynaptic rate. This is given by the weight between the neurons 
times the postsynaptic gain. The gain is the slope of the input-output function (Eq. 2), which is 

monotonically increasing with the postsynaptic cell's firing rate. That is, the steady-state equation 

.1 1 ^—^ 

is Vi = a{J2j JijVj + 9iT, so = na{J2j JijVj + 9iT~^Jij = nany. " J^j. For small a, yi « agf 

is small, and hence the gain is small. In this regime the network is essentially feedforward driven, 

with small modifications by the weak effective recurrent connections. Since individual cells respond 

supralinearly to their inputs, the network sums responses supralinearly. 



3.2 Scaling for large a 

For sufficiently large a, the series in Eq. 8 will explode rather than converge. We expect this 
to occTir for a = Oil). Physically, this occurs approximately when the effective weights become 
strong enough that the excitatory subnetwork by itself becomes unstable in the absence of dynamic 
feedback inhibition. Another way to say this is that inputs are raised to the power n > 1 to 
produce responses which feed back in as inputs; once inputs are sufficiently large, this process is 
explosive, like a nuclear reaction going critical. If the network nonetheless remains stable, it must 
be dynamically stabilized by feedback inhibition (Ozeki et al. 2009, Tsodyks et al. 1997). 

To be dynamically stabilized, the dependence of Q!(Jy-|-g)" on the leading a must be cancelled, 
because otherwise y ~ a, which enters into Jy and is raised to the n*^ power to give^ y ~ a"""*"^, 
which enters into Jy, and so on - the infinite series in powers of a results, which will blow up 
for sufficiently large a. To cancel the leading a, it must be the case that, to leading order in a, 
Jy + g ~ Q:~"- This in turn requires that, to leading order, y has the same a-dependence as g, 
y ~ a'' , so that the leading order of y can cancel the g term leaving only terms of order a" n . We 
define 

/3 = a-n (10) 
Writing y = yo + /3yi, we find the requirement Jy -|- g ~ /3 yields: 

yo = -J~'g (11) 

-J"'g + /3yi = (Jyi)-" (12) 

^This can be done by expanding each power in Eq. 8 as a(g -|- aJ(. . = a(g " + ng''"~^^ * qJ(. . " -|- 

r!i(!^zi)g ("-2) ^ (aJ(. ..)■"). ^ -|- . . ., where .* indicates element-by-element multiplication of two vectors to create 

another vector, and then collecting together the terms of each given order in a. 

^This could be avoided if Jy = to leading order in a, but that requires fine tuning, i.e. it requires Det J = 0. 
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The latter equation shows that yi itself has some further dependence on p. 

These arguments can be translated in terms of r. Once c is sufficiently large, stability requires 
cancellation of the linear dependence of Wr + eg on c, because otherwise r ~ c" which enters 
back into Wr and is raised to the n to yield c" dependence, and so on. Cancellation requires 
that, to leading order, r ~ c, which in turn requires that to leading order Wr + eg ~ C". Writing 
r = ero + e^ri, we find that tq = Ayo = -W~^g and ri = f ^yi = — ^yi, with ri satisfying 

-W-ig + e-'^ri = (Wn)". 

These solutions show that, if the network dynamically stabilizes, its responses arc a sTim of terms 
that are linear and sublincar in the feedforward inputs, that is, responses can add sublinearly. In 
studies of 2-dimensional systems (one excitatory and one inhibitory population), we will find that, 
when the excitatory-neuron element of — J'~''^g is negative, the sublinear term becomes dominant (as 
it must: (yo)^; < 0, so one must have P{y\)E > |(yo)E| for yE > 0) and network behavior becomes 
strongly sublinear. In this case, excitatory firing rates eventually peak and then are ultimately 
pushed to zero with increasing e, i.e. with decreasing j3, but there is a large dynamic range of e 
beyond the supralinear-to-sublinear transition before this peak occurs (see Figure 2). The behavior 
from e = until somewhat beyond the peak yields behavior much like that seen in biology, and so 
we guess that this dynamic range represents the dynamic range of the feedforward input to cortex. 
This will be discussed in Sections 5.2-5.3. 

A more systematic account of the large a (small /3) case can be obtained by formulating a 
solution like Eq. 8 for small a. When the elements of y are > (more precisely: when the elements 
of Jy -I- g are > 0), we can rearrange Eq. 7 for the steady state as 

y = -J-ig + /3J-V-" (13) 
Then we can formally write a steady-state solution as 

y = -J"'g + /3J-'(-J-'g + /3J-'(. . .)-^)-^ (14) 



r = -eW-ig + — W-i(-eW-ig + ^W-i(. . .)-^)-^ (15) 

If quantities in parentheses are positive, a series solution in powers of /? can be obtained from 
Eq. 14 in the same manner as outlined for Eq. 8. When this series converges, which it will for small 
enough (5, it gives a steady-state solution. However, if elements of J~^g are negative, then for small 
enough /? the elements in parentheses will no longer be positive (and correspondingly, as mentioned 
above, in the 2-D case ije is pushed to zero with decreasing j3 at finite /3, so that Eq. 13 fails at 
that point). We can instead regard Eq. 13 as an iterative scheme, y\p+ 1] = — J~^g-|- /3J~^y[p]", 
beginning from some initial condition y[0] (Eq. 8 can also be regarded in this way). Writing this 
as y[p + 1] = /(y [p]), if all of the eigenvalues of the Jacobian of / at the fixed point have absolute 
values less than 1, then the iteration will converge to the fixed point within some basin of attraction 
about the fixed point. Hence with suitable initial conditions, one can find solutions through this 
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iterative scheme, although not for less than that at which some elements of y are pushed to 
zero. 

These scaling arguments provide key insights into the supralinear network (Eq. 3) that is con- 
firmed by other analysis and simulations: for small a, recurrence is weak and the network supra- 
linearly adds responses to different feedforward inputs; with increasing a, there is a transition, for 
a 0(1), to a dynamic stabilization that leads responses to add sublinearly. Note that individual 
neurons still supralinearly sum the net (feedforward plus recurrent) inputs they receive, but the 
network "conspires" to deliver net input that is so strongly sublinear that, even after the neuron 
raises its net input to the power n, its responses add sublinearly. We have found in both high- 
dimensional and 2-dimensional simulations, and we will show below for the 2-dimensional case, that 
stabilization will occur provided feedback inhibition is sufficiently strong and the inhibitory time 
constant is not too slow relative to the excitatory time constant. This transition from supralinear 
to sublinear behavior in turn appears to undcrly a wide variety of nonlinearities in neocortical 
behavior. 

3.3 Comparison to the balanced network 

Van Vreeswijk and Sompolinsky (1996, 1998) introduced the "balanced network" model (see also 
Renart et al. 2010). They considered a circuit of stochastic excitatory and inhibitory units that 
could have firing rate or 1, and studied the conditions in which the network would dynamically 
find its way to a balanced state in which the mean input is subthreshold, yet firing rates are nonzero, 
meaning firing is driven by fluctuations. They assumed each unit received K inputs of strength 

or a net input of strength ^/K, for K large {e.g., thousands of inputs). The mean field equations for 
the average E and / firing rates are the 2-dinicnsional version of the rate equation, Eq. 1, for one 
E and one I population, where both W and eg are of order y/K,^ and the function / is a sigmoidal 
function rising from to 1 as the input moves from approximately —3 to 3, and saturating at or 
1 for smaller or larger values respectively. To be in the balanced state, the mean firing rate must 
be neither nor saturated at 1, so the net input must be 0(1) (i.e., between —3 and 3). 

Thus, the condition for the balanced state is that Wr + eg ~ 0(1) where both W and eg are 
0{\/K). The solution, much as in our scaling argument, is to write r = fq -|- :^ri -|- . . ., where tq 

and ri are 0(1) and the dots represent higher-order terms in The balanced condition is that 

the 0{-\/K) term in the input vanishes, that is, Wro = —eg, leaving as input only the 0(1) term 
^^ri and terms that are 0{^=). 

^/ K y K 

The condition for dynamic stabilization to achieve the balanced state, fq = — eW~^g, is of 
course identical to the condition we have found for dynamic stabilization in our model. ^ Although 
the condition is formally identical, the meaning is different in crucial ways: 

■^The input is expressed in units of a variance term which itself is dynamically determined, but this term is 0(1) 

and docs not impact the points made here. 

^We wrote the leading term as cro rather than ro, but the leading terms are identical. 
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1. In the balanced model, the stabilization is required because inputs are large and must cancel 
to leave something small. In our model, the stabilization arises because the dynamics are 
explosive, due to the supralincar input-output function, and occurs when none of the inputs 
arc large. This can be seen by recalling that in our model the effective recurrent weights are 

knW and the effective input strength is k^c, and that a = (^fc"||W||^ (ji^cj . In the 

balanced network, the recurrent weights fe»||W|| and the feedforward weight k^c are both 
O(v^), so a is 0(^2 ). In contrast, stabilization arises in our model when a is 0(1). 

2. In the balanced network, the second-order term is negligibly small relative to the first- 
order term rg (because the stabilization is to cancel large things, leaving something small). 
The first-order term is linear in the input, tq = — cW~^g, and so in the balanced network 
responses are always linear in the input. In our model, because the stabilization occurs when 
inputs are small, the second-order term can be comparable to or larger than the first-order 
term over a wide dynamic range, enabling a variety of sublincar behavior. 

In particular, in our model elements of tq can be negative, meaning that for such an element 
ri > |ro| over the relevant dynamic range of behavior (discussed in more detail for the 2-D 
model in Section 5.2). In the balanced model, since all terms except tq are negligible, the 
elements of fq must be positive for activity to be nonzero. 

In sum, in the balanced network, inputs are huge relative to the distance from rest to threshold, 
and must dynamically cancel for the network to neither saturate nor have activity but instead be 
in the fiuctuation-driven regime. The dynamic cancellation or stabilization yields network responses 
that are always linear in the input. In our model, the supralinear input-output function renders 
the network explosive - input is raised to a power greater than 1 to produce responses, which feed 
back as input. Stabilization against this explosive nonlinearity arises when inputs are relatively 
small, yielding a range of sublinear behavior. 

4 Reduction to a 2-dimensional system 

Most of our analysis hereafter will focus on a 2-dimcnsional system of one excitatory and one 
inhibitory population, as it is difficult to say much in general in higher dimensions. A 2-D system of 
one E and one I population can be derived as a mean field equation from higher-dimensional models 
in which E and I neurons have random connectivity {e.g. Renart et al. 2010, van Vreeswijk and 
Sompolinsky 1998). In particular, if the high-dimensional model involves integrate-and-fire neurons, 
their input-output functions in the fluctuation-driven regime can be reasonably approximated by 
power-law functions (Hansel and van Vreeswijk 2002). 

Here we consider a higher-dimensional system with structured connectivity. We show a heuris- 
tic derivation of a 2-D system that preserves a surprising amount of the behavior of the higher- 
dimensional system. We then show how the conditions for "normalization" - sublinear addition of 
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responses to multiple stimuli - in the high-D system can be expressed as simple conditions in the 
2-D system on the growth of r with increasing c, or the growth of y with increasing a. 



4.1 Reduction 

We consider a topographic network, with pairs of excitatory (E) and inhibitory (I) units arranged 
on a 1-D or 2-D grid with periodic boundary conditions. Stimuli are localized on the grid, though 
there may be more than one localized stimulus present. For simplicity, we assume a single time 
constant for all E cells and one for all I cells. We let represent the position on the grid, and rE{0) 
and rjiO) the excitatory and inhibitory firing rates at position 9. Thus, we can write Eq. 3 as 

TE^^I^ = -rE{e) + k{i^JEE*rE{e)-^I^JEi^ri{e) + cgE{e)y'' (16) 

r/^^ = -ri{e) + k{^JiE*rE{e)-^Jn*ri{e) + cgi{e)y'' (17) 

Here, Jxy * ry(e) = JxY{0,9')rY{9'). 

We will consider "normalization" , the sublinear addition of the responses to two stimuli (Caran- 
dini and Heeger 2011). We let one stimulus be centered at 9 = 0. We let jxY be the vector of weights 
JxviO, to position 0, normalized so that the weight from position is 1, and let Jxy = ^xy(0, 0), 
the actual value of the weight from position 0. Similarly, we let te-, r/ be the vectors of excitatory 
and inhibitory firing rates, respectively, normalized to equal 1 at position 0. Then the equations 
for the units at position are 

^^drs^ = -rB(0) + k (^^JEETEiO) (jee ■ ve) - ^JEiriiO) (jei ■ r/) + c5b(0)) " (18) 

r/^^ = -r/(0) + k (^pJiEVEiO) (jiE ■ te) - ^Jiiri{0) (jn ■ rj) + 057(0)) " (19) 

Although we had previously incorporated changes in |g| into c, we now take addition of a second 
stimulus to simply alter g, so that in particular addition of a second stimulus that gives no input 
to position does not alter ^^(O) or gi{0). 

We now make the ansatz that, as the stimulus changes, the four dot products are all scaled by a 
common factor, whether the stimulus changes in strength (changing c) or in shape (changing g, e.g. 

adding a second stimulus). We can then write Jxy = Jxy ^r^^ and ^ = ■;/; (Jee • t^bV where ^' 

[jEE-rE) \ J 

includes the common scaling factor (note that the matrix J of these J's need not satisfy || J|| = 1). 
Since the weights jxY are non-negative, the effect of adding a second non-negative stimulus is to 
increase ^, hence "normalization" of the E or I population corresponds to a decrease in firing rates 
rE{0) or r/(0), respectively, with increasing '^^^^^ < 0, '^^^^^ < 0. 

With this ansatz, we obtain 2-dimensional equations. Letting ve = '"£;(0), gE = ^^(O), etc., we 
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obtain the equations 



drE 



—TE + k {^Jeete — ^Jeiti + cgE) 



.n 



(20) 




ri + k {"^JiETE - '^Jnri + cgi) 



.n 



(21) 



This is just the 2- dimensional version of Eq. 3, and is equivalent to Eq. 6 for 2-dimensional y with 
3' replacing il) in the definitions of y and a (Eq. 5). 

The ansatz, of course, is not in general true, but it can be close enough to true to give a 
good account of the higher-dimensional system. To illustrate this, we simulate the model on a 
one-dimcnsional ring, which we can think of as representing preferred orientation.^ We consider 
180 E/I pairs at grid positions separated by 1° in preferred orientation, with 0° = 180°. All four 
connection types have the same width, following evidence that excitatory and inhibitory inputs 
received by cells in upper layers have similar orientation tuning (Anderson et al. 2000a, Ferster 
1986, Marino et al. 2005, Martinez et al. 2002). The connectivity takes the form 



where d{6,6') is the shortest distance around the circle between 6 and 9', Jee = 0.0441, Jje = 
0.04158, Jei = 0.0231, and J// = 0.01827. (Here we arc not sticking to our convention of ||J|| = 1, 
or in other words, we are taking ip = 1 with these values for J.) We take Cori = 32°, but also 
consider other values in Fig. IC. Other parameters are te = 20 ms, tj = 10 ms, k = 0.04, and 
n = 2.0. We consider stimulation by either one oriented luminance grating or two orthogonal 
gratings. Each grating is represented by a Gaussian-shaped curve of feedforward input with width 
(standard deviation of the Gaussian) 30° and height c; a single grating is centered at = 0°, a 
second added grating is centered at 9 = 90°. We also consider varying the width of the feedforward 
Gaussian, for a grating centered at 9 = 0°. For any given stimulus (1 or 2 stimuli, stimulus height 
c, given stimulus width) the equivalent 2-D model is found as follows: we use the same J's and 
other parameters, and take ^ to be the value of the convolution, at 9 = 0, of the connectivity 
Gaussian (Eq. 22 with Jxy = 1) with the input pattern for c = 1, which we use as a surrogate for 
the shape of the response. For the default Gaussian width, this yields * = 54.86 for one stimulus 
and ^ = 67.67 for two stimuli. 

The result is that the reduced 2-D model accurately reproduces the behavior of the full model 
(Fig. 1). The firing rates of the cells at ^ = vs. stimulus strength closely match the firing rates 
in the 2-D model (Fig. lA). Both models show a similar transition from supralinear summation 
of responses to the two gratings for weak stimuli to sublinear summation or "normalization" for 

®The paradigm we study here supjjrossioii of response to one orientation by presentation of an orthogonal 
orientation - is known as "cross-orientation suppression". In VI, this appears to be primarily mediated by sublineax 
addition of the feedforward inputs to VI evoked by the two stimuli (Lauritzen et al. 2001, Li et al. 2006, Priebe and 
Ferster 2006). However wc use this paradigm to study how the model cortex sums responses to multiple stimuli, 
assuming the feedforward inputs sum linearly. 



WxYi9,9') = JxYe 



(22) 
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Figure 1: 

Two neuron approximation (right) of the full ring model (left). (A) The reduced version of tlie 
model produces qualitatively siruilar curves of response vs. stimulus strength c (for the full model, this 
is the response of the cells at 6* — 0). (B) Full and reduced models show a similar stimulus-strength- 
dependent transition from supralinear summation (weight > 1) to sublinear summation (weight < 1) of 
the responses to two gratings. For the full model, for either E or I cells, we let Ri{9), R2{d), and Ri2{d) 
be the response to one grating, the other grating, or the superposition of the two, and find the weight 
w that gives the least-squares-fit to Ri2{d) = w {Ri{0) + R2{d))- For the reduced model, the correspond- 
ing responses are i?2, and i?i2, and the weight is w = ^^^^ ■ (C) Full and reduced models have 
nearly identical stimulus-strength-dependent tuning for the width of a feedforward stimulus (full model: 
width of Gaussian stimulus centered at 9 = with given stimulus strength c that gives strongest re- 
sponse in cells at = 0; reduced model: '5 is computed for each stimulus width, and plot shows width 
whose ^' gives maximal response). In all curves, red shows E cells and blue shows I cells. All responses 
are steady-state responses. Full model solutions found by simulating until convergence to steady state. 



stronger stimuli (Fig. IB). The network also shows a form of surround suppression, in which the 
"summation field size" - the stimulus width that yields maximal response for a given stimulus 
strength - shrinks monotonically with increasing stimulus strength, as is well known in real space 
(rather than orientation space) for VI cells (Cavanaugh et al. 2002, Sceniak et al. 1999), and this 
behavior is extremely similar in the full and reduced models (Fig. IC).^ Thus, the 2D model can 
provide a good basis for understanding more general models. 

4.2 Conditions for Normalization 

We consider steady-state r or y and use expressions like to refer to the dependence of the steady 
state on parameters. We have seen that the condition that rx {X G {E, I}) exhibits normalization 
in the high-D model is equivalent to the condition < in the 2-D model. Since Eqs. 20-21 are 
equivalent to Eqs. 5-6 with * replacing ip, we revert to the notation of Eqs. 3-6 and use ■0. 

We work with the 2-D model and express the conditions < as a single vector condition. 

We note first that $ = c'^ = c (ig - |i) and J = = kc^'''^. Putting these together 

we find ^ = ^ (a^ ~ a) • Thus, the condition for normalization is that y grow more slowly than 

linearly with increasing a: ^ < ^ or < 1 or, roughly, that y '-^ for p < 1. As we have seen, 
p becomes less than 1 precisely when the transition from the supralinear to the sublinear scaling 
regime occurs. 

We can reexpress this in terms of r. Using algebra similar to the above, we find ^ = 
+ J^)^ from wliich we find that ^ < ^ is equivalent to | < n^. Thus, the 
condition for normalization is that r grow more slowly than with increasing c: ^|^^ < n or, 

''Note that this "summation field size" for orientation selectivity should not be confused with the orientation 
tuning width, which is the width of the orientation tuning curve obtained by studying response vs. single orientations 
(more precisely: studying response vs. center orientation, using stimuli that evoke a fixed curve of feedforward input 
vs. orientation that is symmetric about the center orientation). The orientation tuning curve, representing the set of 
single orientations that can drive the cell, is analogous to the "minimal response field" in real space, which represents 
the sum of the set of small regions in visual space in which appropriate light stimuli can evoke spiking responses. 
The minimal response field in real space is invariant with stimulus contrast (Song and Li 2008), and so too is the 
shape of the orientation tuning curve (Anderson et al. 2000b, Ferster and Miller 2000, Skottun et al. 1987) (contrast is 
monotonically related to the firing rate of the inputs to cortex {e.g. Ohzawa et al. 1985)). The fact that the summation 
field size in real space is larger than the minimal response field indicates that stimuli in regions where light cannot 
directly drive spikes can facilitate responses to stimuli in the minimal response field. Recall that the size of this 
facilitating area shrinks with contrast. The model suggests that the same may be true in the orientation domain, 
in terms of cortical processing of feedforward input to cells of different preferred orientations. However, attempts to 
test this idea will likely be compromised by two facts: (1) simultaneous presentation of multiple orientations does not 
yield linear summation of the input to cortex evoked by the individual orientations (Lauritzen ct al. 2001, Li ct al. 
2006, Priebe and Ferster 2006) and (2) varying the feedforward orientation tuning by changing stimulus attributes 
- e.g. a sinusoidal luminance grating of a given size provides drive to cortical cells with an orientation tuning that 
narrows with increasing spatial frequency, and similarly a longer bar drives narrower orientation tuning than a shorter 
bar - also changes other attributes to which the neurons are independently sensitive, such as spatial frequency or bar 
length. 
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roughly, that r for p < n. Again, p becomes less than n at the transition from supralinear to 

sublincar scaling. 

Finally, noting that the steady state condition is r = fc(Jr + eg)'", without loss of generality we 
write Jr = cf(c) for some vector function f of c, so that the steady state condition r = fc(Jr + cg) "' 
becomes r = c"(f (c) + g) ". Thus we see that a component of r grows more slowly than precisely 
when the corresponding component of f(c) is a decreasing function of c (that is, for corresponding 
components r and f ■, % < precisely when /'(c) < 0). Thus, the condition for normalization can 
alternatively be expressed as the requirement that Jr grow more slowly than linearly with c, i.e. 

dc ^ c dine ^ 

5 Analyses of the 2-Dimensional Network 

We will assume throughout this analysis that qe > 0, gi > 0. We will use the following definitions: 

= BetJ {-3-^g)j^ = JngE- JEigi (23) 
ni = Bet 3 {-J-'^g) J = JiEgE- JEEgi (24) 

We also note that there are three possible conditions: (1) (— J~^g)g > and (— J~-'^g)^ > 0; (2) 
(-J-^g)g < and (-J-^g)^ > 0; and (3) (-J-^g)^ < and (-J-^g)^ < 0. The 4**^ condition, 
(— J^^g)^ > and (— J~^g)^ < 0, is not mathematically possible for gE > and > 0: Q,e > 
and <0 together imply Det J < 0, and similarly fi^; < and ilj > together imply Det J > 0. 



5.1 When Does the Network Dynamically Stabilize? 
5.1.1 The case of infinitely fast inhibition 

We first analyze the case of infinitely fast inhibition, ti/te = 0, with constant feedforward inputs. 
We show that in this case the network, if Det J > 0, the network is always driven to a stable fixed 
point from arbitrary starting conditions. The condition Det J > means that feedback inhibition 
is sufficiently strong: JeiJie > JeeJii- In addition, we show that if Det J < 0, sufficiently large 
initial firing rates will cause the system to "blow up", i.e. firing rates will grow arbitrarily large. 

With r/ = 0, the value of yj is "slaved" to, or instantaneously set by the value of, tje according 
to the ^ part of Eq. 6 for y. Because of the nonlinearity, we cannot solve this for yj as a function 
oi yE, but we can instead solve for as a function of yj: 

Substituting this in the part of Eq. 6 yields, after a bit of algebra, an equation for induced 
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by the slaving of yi to the ue dynamics: 



dyi _ nanyj"" I (yi\h " + t ,7 (yi\^ 



l + Jiinanyj" \ "^lE \ / / 

n ^^^^ 

For sufficiently large yj, if Det J > 0, the term inside the parentheses in the ° ^ (. . .)'" term 

■■'lE 

will be negative, and so will be set to zero after the thresholding involved in the () " operation. 
The dominant term will then be the —Jiiyi term, which is negative. So for sufficiently large yj, 
^ < 0. On the other hand, if Det J > 0, then for sufficiently large yj, the (. . .)'" will be positive 
and larger than the sum of the other terms, so that ^ > and, since increasing yj will increase 
^ > 0, this derivative is ever-increasing. 

For sufficiently small yj, ^ > if either gj or gE is nonzero, which can be seen as follows. 
For sufficiently small yi, the source terms gi and fie, if nonzero, dominate the terms involving y/. 
Both gj and the (. . .)'"" term containing ^e are non-negative, so if either is positive ^ will be 
positive; if Qe > 0, the (. . .)'^ term is positive; if ^Ie < 0, this implies gi > (given that at least 
one of gi and is nonzero, and that both are non- negative) . 

Thus, for Det J > 0, yj is driven to a stable fixed point, and yE is then determined from Eq. 25, 
so the system will arrive at a stable fixed point. Note that the system could have multiple fixed 
points with varying levels of yj. The topology of flow along the yj axis tells us that there must 
be an odd number of fixed points, alternating from stable to unstable to stable with increasing yj, 
with the outermost fixed points (those with lowest and highest yi) being stable. In the simplest 
case, there is a single stable fixed point. In addition, for Det J < 0, the system will blow up for 
sufficiently large initial firing rates. 

5.1.2 More general requirements for stability 

Changes in the time constants can alter the stability of the fixed points, but do not alter the number 
or positions of the fixed points. The results of the previous section tells us that, for Det J > 0, the 
system always has a fixed point that is stable for tj = 0. We consider such a fixed point, and ask 
when it retains or loses stability for finite tj. 

We let the fixed point be ( J , and assess stability by linearizing the dynamics about 



this fixed point. We let q = tj/te > 0. Setting r = in Eq. 6, the matrix T is given by 




Define the matrix ^ = nan [ ^ ^ | . Writing the identity matrix as 1, 



1 

q 

the Jacobian matrix of the 2-D system is: 

J = T-^ (*J - 1) (27) 
A fixed-point of the dynamics will be stable if J has a negative trace and a positive determinant. 
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The negative trace condition is Jee < <Jii, which becomes 



1 



Jee = na^yj/ Jg^ - 1 < OR \ Jee > AND q < ""'"^^i "^^^^^ ) (28) 

nany^"" Jee - 1, 

The condition Jee < means that the excitatory subnetwork by itself is stable (or marginally 

stable), which guarantees that the network will always be stable, since only excitatory instability 
can destabilize the network. When the excitatory subnetwork is unstable, we can further reduce 
the condition on g for n = 2:^ 



Jl + 4:aJn {gi + JiEyE) , „ ^, , ^ 

9 < TTj ^= — ('^ = 2, ^\/ay^JEE > 1) (29) 

^JEE^/ayE - J- 

The determinant condition, Dei J > 0, is always true for any fixed point that is stable at 
q = 0. To see this, note that the sign of the determinant does not depend on for g > (because 
Det AB = Det A Det B for any matrices A, B, and Det — q)- So if we prove that Det J" > for 
arbitrarily small > 0, we will have shown that it holds for all q > 0. For q = 0, the determinant, 
which is the product of the two eigenvalues, was infinite: because the fixed point was stable, both 
eigenvalues had negative real part: one real part was infinite, corresponding to the instantaneous 
flow onto the inhibitory nullcline (the line in the yE/yi plane on which ^ = 0); the other was 
finite, corresponding to the flow along the nullcline converging onto the fixed point. (Since the 
two real parts were unequal, both eigenvalues were real.) As q is moved inflnitesimally from 0, 
the inflnite eigenvalue becomes a large but finite negative eigenvalue, while the finite eigenvalue is 
perturbed by arbitrarily small amounts as q is made arbitrarily small. This means that there is 
a range of g > for which the eigenvalues continue to have negative real parts, and therefore for 
which the determinant condition holds. Therefore, the determinant condition holds for all q. Thus, 
for a fixed point that is stable for g = 0, the fixed point remains stable so long as condition 28, or 
condition 29 for n = 2, is satisfied. 

We also note that, for the case n = 2 and for g < 1, a sufficient condition to conclude that 
there is only a single fixed point, which is stable, is Det J > and Jj^^ < JieJii, which can be 
seen as follows. The determinant condition is Det (*J — 1) > 0. We note that, for an arbitrary 
2-dimensional matrix M, Det (M — 1) = DetM — TrM + 1. Thus, the determinant condition is 
Det $J > Tr $J — 1. Since Det $ > (because firing rates and a arc > 0), this condition will be 
satisfied if Det J > and Tr*J < 1. The trace condition for stability is TrT~^$J < 1 + g. But, 
for g < 1 and given the structures of * and J, TrT~^*J < Tr*J, so the condition Tr*J < 1 
ensures that the trace condition is also satisfied. This condition is 

n — 1 n — 1 

JEEVif' - JiiyT^ < — r (30) 

nan 

®This condition is found by solving y/yj = y/a{JiEyE — Jiiyi + Qi) as a quadratic equation for y/yi. Discarding 
the negative solution, this yields = ^+V^+'^^^^'^'^+-^ieve) ^ Substituting this into Eq. 28 for n = 2 yields 
Eq. 29. 
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For n = 2, we substitute the solution for ^fyj as a function of yE (footnote 8) into Eq. 30 for n = 2 
to find Jgg — J IE J II < ^^tays" ' right side is positive, a sufficient condition for this to 

be true is < JieJii- Recall that, if there is more than one fixed point, some will be unstable 
at g = 0, and they must remain unstable for some region of small but finite q. Since this condition 
guarantees that any fixed point is stable, we conclude that there can only be one fixed point, which 
is stable, when this condition holds. 

In summary, for q = ti /te = 0, the network always flows to a stable fixed point if Det J > 0. 
For q > 0, a fixed point that is stable at g = remains stable when Eq. 28 or, for n = 2, Eq. 29 
is satisfied. Note that this condition does not ensure that the network always fiows to a stable 
fixed point; for nonzero q there may be initial conditions outside the basin of attraction of the 
stable fixed point or points. A condition that ensures that any fixed point is stable for q < 1, and 
therefore that there is only one fixed point, is Det J > and j'^^ < JieJii- If this fixed point is 
the only attractor (there are no limit cycles), this ensures that the network will flow to the stable 
fixed point. Excepting Eqs. 28-29, these conditions involve feedback inhibition being sufficiently 

strong: JieJei > JeeJii, and Jie > -f^- 

In Fig. 2, bottom row, we will illustrate the range of g's yielding stability for various parameter 
choices with n = 2. 

5.2 The case (— J~^g)g < and supersaturation 

We consider Eq. 3 for r, but substituting V'J for W. We restrict to the case Det J > 0, which ensures 
a stable fixed point for at least some range of ^ > 0. We note that for Det J > 0, (— J~^g)^ < 
and (— J^-'^g)^ < arc equivalent to VIe < and flj < 0, respectively. 

We shall equate increasing or decreasing c with increasing or decreasing stimulus contrast. This 
is based on the fact that the contrast of a visual stimulus is monotonically (but nonlinearly) related 
to the firing rate of the inputs to VI from the lateral geniculate nucleus (LGN) {e.g. Ohzawa et al. 
1985). 

In simulations, we find that if il^; < 17/ < for qe = gi, then te grows with c for a range 
of c considerably beyond the transition from supralinear to sublinear behavior, but ultimately 
peaks and is pushed back to with increasing c (see Fig. 2A). The inputs to cortex have limited 
dynamic range {e.g., stimulus contrast cannot increase beyond 100%), and so we imagine that 
this circuit may model cortex but that the maximal input strength seen biologically cannot drive 
excitatory responses too far beyond their peak. The decrease in response with increasing contrast 
after a peak response is referred to as "supersaturation", and is seen in virtually all VI cells for 
contrasts larger than about 75% (Ledgeway et al. 2005, Li and Creutzfeldt 1984, Peirce 2007, Tyler 
and Apkarian 1985). This model behavior provides one possible explanation for supersaturation, 
although supersaturation might also in part refiect a supersaturation of inputs, e.g. if feedforward 
inhibition (Bruno 2011) overtakes feedforward excitation with increasing contrast. 

Here we analyze this behavior. We shall find that (1) if r is a stable fixed point, then ^ and 
^ are negative precisely when ^e < and 0/ < ^' „_i , respectively (and so in 
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particular can only be negative if < or $7/ < 0, respectively); (2) if fi^; < 0, then there is a 
stable fixed point with = at a finite positive value of c, which we calculate; and, (3) for n = 2, 

2 

if in addition Ue < ^^i, then the set of fixed points rE{c) reaches a maximum for increasing c 

with ^ = before being pushed to zero, and we calculate the corresponding c and peak value of 
rE- 

The condition 0,e < states that the linear term in c in the high-c expansion for (Eq. 15) 
is negative, driving to zero; while the requirement in condition three states that the linear term 
in the expansion for r/ is either positive, or not so negative as to disrupt the ability of inhibition 
to drive to zero. 

These results suggest, but do not prove, that will be driven to zero for arbitrary n whenever 
CIe < (although there is a stable fixed point with rg = at finite c, we have not proven that it 
is the only stable fixed point). We note that r/ can never be zero for finite c if 57 > or r^; > 0, 
so even for < 0, rj can never be driven to zero with increasing c. 

For Qe < < 0, gE = gi, we find in simulations that r/ only increases with increasing c, 
and that the values of c at which ve goes to zero and, for n = 2, at which peaks and the 
corresponding peak all are as calculated (see Fig. 2A). We speculate that, fov CIe < and 



< 0, while rr always becomes large enough to set < and so ultimately to drive to zero. 

When Cli < Qe < for gE = gi, we find unbiological behavior in simulations in which both 
rE and r/ jump to very high levels at very low c, after which rE monotonically decreases and is 
ultimately pushed to (see Fig. 2E). Numerical calculations suggest a discontinuity at the jump, 
which may explain why our calculations do not find a zero of ^ for real positive c in this case. 
We have not tried to analyze this behavior. 

5.2.1 When can or rj decrease w^ith contrast? 



Ue < 




^i, where f{n) may equal n or may equal 2, rg can never become large enough to set 




V rj" J 
(1 — ^#rJ)~ *rg, which gives 



Then a simple calculation shows that ^ 



dc 



*r(V'J|+g) or 



dr ^ 
dc 



dr 

dc 




(31) 



Det (1 - V*rJ) 
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Stability requires that Det (1 — *rW) > 0. Thus, this expression shows that, for a stable fixed 
point, ve oi ri decrease with contrast precisely when 



< (^<0) (32) 

ntpk " ^ ^ 

5.2.2 The c at which te becomes 

The c > at which te first becomes with increasing c can be determined as follows. First, at this 

c, rj = cgE/ipJEi, because this is the value of rj that sets the input to to zero when = 0. 

The equation for the r/ steady state then yields ^j^^ = k [cgi — cgEj^^ = kd^ ( ^gf ) ' '^^^ 
right side gives zero unless rj^; < 0, so a solution for c 7^ exists only for fi^ < 0. In this case, one 

can solve to find c = Jei ( u u^r, \n \ ■ This corresponds to a = f or i3 = ^ — j-.^ 

Note that any fixed point ue = 0, yj > is stable for any q since the Jacobian matrix is 

i7 = na^ 1 1 / \ which has two negative eigenvalues (equal to 

the two diagonal entries of J'). 

This shows that ve = 0, rj = cgE/ipJEi is a stable fixed point for this value of c, but does not 
rule out the existence of other fixed points. 



5.2.3 Peak firing rate and corresponding contrast 

Prom the steady-state equation for r^;, we can solve for rj in terms of r^;. Substituting this into 
the steady-state equation for rj, we obtain an implicit equation for the steady-state firing rate of 

TE- 

-A 



1 

TE 



V'Det J \ \ k J \ ipJEik 



(34) 



We now restrict to the case n = 2. We take the derivative of both sides w.r.t. c, yielding an 
implicit equation for ^ which we solve for ^ . We then solve for ^ = to find the c at which 



®Once rE has been pushed to zero, for increasing c, rj continues to increase according to ri = k{cgi — ipJiiVi)^ , 
which for n = 2 has the solution n = ^ ikiii'^j^ 1 ^"^^ remains 0. 
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a maximum firing rate can occur: 

Jee^te 




If we then substitute c^^^^ for c in Eq. 34, we can solve explicitly for the corresponding maximum 
value of rg, r^*^. For this to have a real value - that is, for a maximum of te as a function of 
c to exist - it must be the case that (— J~^g)g < 0. We restrict to the case Det J > 0, which is 
required for stability, so that the requirement for a real solution is < 0. We then obtain the 
following expression for r^*^. With $7^; < 0, define: 



glni + 2g'j\^E\-'^gi^\^E\ {glfli + g]\nE\) 

(36) 



Then the maximum firing rate r^"^ is^^: 



(38) 



Note that for this to be real, an additional condition needs to be satisfied, namely g^^i > ~9'j\^E\- 
Substituting this expression back into Eq. 35, we obtain an expression for c^^^ in terms of 
model parameters: 

Explicit calculation shows that this is always a maximum: the 2"*^ derivative < 0.^^ is 
not guaranteed to be positive - this is governed by rather complicated conditions on the g^s and 
J's - but in practice we have found it to be positive for the simulation parameters we have used. 

When c^^^ is positive, then reaches a maximum as a function of c at c™*^^. Equations 37 
and 38 show that, in this case, both the maximum excitatory firing rate that can be achieved 
by the network and the contrast at which this maximum is achieved decrease with increasing tp. 
Thus, when the 2-D reduced model, Eqs. 20-21, accurately captures the high-dimensional model, 
Eqs. 18-19, then in the high-D model, if the stimulus is widened or a second stimulus is added, the 
maximum excitatory firing rate will go down and will occur at a lower contrast. 



These calculations were done in Mathematica. We solved for the numerator of the expression for ^ being 0, 
assuming that the denominator was not simultaneously zero. We calculate further that when the numerator is zero 
(c — c™"^). and for Qe < 0, the denominator is zero when te ~ -r, — ■>, , ^^'^^'^ — ^ ^ . So long as this does not 
coincide with the value of r_B for c — c™^^ (Eq. 37) the procedure is fine. 

^^Thcrc is a second solution in which the sign in front of the square-root term in k (Eq. 36) is positive. For 
parameters we have used in simulations the first solution gives positive c and the second gives negative c, so we have 
focused on the first solution, but there may be parameters for which this situation is reversed. 

^^For this 2nd solution mentioned in footnote 11, the 2nd derivative is positive precisely when Clj > 0. 
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^max correspond toa=^ ('M^ - Jeek + 2y^) and ve = J^r^"^ = y-^^ ^ y 

However, note that this is not a maximum of the i/e vs. a curve, but rather occurs for a higher than 
that peak, where the curve has a negative slope. We saw in section 4.2 that ^ = + j > 

so ^ = implies ^ = — (^^rryy^) i-^- Y is locahy evolving as a~^'^~^\ 

In sum, for Det J > and n = 2, the steady state solution for rE has a maximum value as a 

2 

function of c, given by Eq. 37, precisely when (1) < and (2) < ^fi/. 
5.3 Steady-state solutions for different petrameter regimes 

In Fig. 2 we illustrate model behavior, as a function of stimulus strength c, for 5 parameter regimes, 
with Det J > 0, n = 2, and qe = gi in all cases. The 5 parameter regimes arc: Jlf; < and 0/ < 0, 
with either 9.e < (Fig. 2A) or Qe > (Fig. 2E); ^e < and Qj > (Fig. 2B); and Qe > 
and Qi > 0, with either Qe < (Fig. 2C) or rig > 0/ (Fig. 2D). We chose parameters relatively 
arbitrarily, by starting with a set of parameters that had worked well in simulations of the ring 
model (Column A) and changing small sets of parameters to change the regime. However in small 
amounts of studies of other parameters in the different regimes we have found behaviors to be 
similar, with one exception. For the case < and 0/ > (column B), the excitatory firing rate 
can peak and be driven to zero without ever reaching a level at which the excitatory subnetwork is 
unstable, as manifest in the figures as stability for all possible values of q (bottom row, described 
below). This occurs because $7/ > is compatible with Jee = 0, and so for weak enough Jee the 
behavior can be similar to Jee = behavior. (For Qe > and 0/ > 0, ve ultimately grows linearly 
with increasing c for large enough c, so given any small but finite Jee the excitatory subnetwork 
must ultimately reach instability.) 

Wc illustrate behavior across a large range of c, sufficient to include the point at which te 
is pushed to zero for cases with Qe < 0. However, we imagine the dynamic range of cortex, 
corresponding to the dynamic range of the firing rates of the inputs to cortex, corresponds to a 
smaller range, e.g. through c = 100 (as in Fig. lA, reduced model, 1 stimulus, which uses essentially 
the same parameters as Fig. 2A). 

For each set of parameters, we first illustrate firing rates (top row), with red and blue indicating 
rE and rj respectively. As expected, parameters with Qe < (columns A,B,E) all show ve 
eventually pushed to zero with increasing c, while those with Qe > (columns C,D) show ve 
moving toward linear growth with increasing c. The combination Qj < Qe < (column E) leads, 
as mentioned previously, to unbiological behavior in which both E and I rates abruptly jump 
(discontinuously, in numerical calculations with c discretized in 0.00001 steps) to high rates at low 
c, after which monotonically falls with increasing c. 

We next illustrate normalization weights (second row), computed just as in Fig. 1, so that 
weights > 1 indicate supralinear summation in the corresponding ring model and weights < 1 
indicate sublinear summation. All but the case Qe > Qi > show a regime of supralinear 
summation for very low contrasts (behavior in all cases is sublinear for c > 10), although the 
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Figure 2: 

Behavior of the 2D Model in Different Parameter Regimes. Each column indicates a different 
parameter set, as indicated. In all cases, Det J > 0, n = 2 and Qe = 9i = The first column uses the 
same parameters as the 2-D reduced model in Fig. 1, except that weights here were rounded to 2 significant 
figures. In all figures the horizontal axis is stimulus strength c. Top row: E (red) and I (blue) firing 
rates, and r/, at fixed point. For cases with $7^; < 0, dashed vertical lines indicate analytic prediction 
for c at which peaks (Eq. 38) and at which te goes to zero (Section 5.2.2). Second Row: Weights 
reflecting supralinear summation (weight > 1) or sublinear summation (weight < 1) in an equivalent ring 
model, computed as in Fig. IB. Again, red and blue indicate E- and I-cells, respectively. Inset in column E 
shows supralinear responses at low values of c. Third Row: Iterative solutions for rE in the low-contrast 
regime. We use (Eq. 39) for ue^ and graph rE[t] — yE[t]c/i' vs. c, where t is number of iterations. Black 
curves are exact solutions; red to yellow represent iterative solutions for 1, 5, 10, 14, 19 iterations. Iterative 
solutions are shown only over the range for which they are real. Fourth Row: Iterative solutions for in 
the high-contrast regime, using Eq. 40. Conventions as in 3rd row, except now blue to cyan represent 1, 5, 10, 
14, 19 iterations. Fifth Row: Values oi q = tj /te separating regions in which fixed point is stable (below 
red line) vs. unstable (above red line). Parameters used: W = -^J, where -0 = 54.86, the value of ^E* in Fig. 1 
for the default Gaussian width. The elements of J are: Jee — 0.044 except in (D), where Jee = 0.014; 
Jei = 0.023; (A) Jje = 0.042; Jji = 0.018; (B) Jie = 0.082; J/j = 0.018; (C) Jie = 0.082; Jn = 0.038; 
(D) Jie = 0.062; J// = 0.088; (E) Jie = 0.039; J// = 0.018. We convert from a to c using k = 0.04. This 
yields the following values of a for columns A to E, respectively: a — .148c, .213c, .226c, .243c, and .144c. 



supralinear behavior is weak for 0/ > 0. 

The third and fourth rows of Fig. 2 ilhistratc the iterative solutions that stem from the scaling 
solutions in the low- and high-contrast regimes. The values of J used (listed in legend of Fig. 2) are 
not normalized to have ||J|| = 1, so for these iterations we take J = J/||J|| and ip = V'||J|| where 
II J|| is the 2-norm of J (the maximum singular value of J), with a = kc^~^'ip. We then show results 
as te = yEc/'4' vs. c. 

The small-a (low contrast) iterations arc shown in the third row of Fig. 2. Here, we treat the 
equation y = a(Jy + g)'*^ as the recurrence relation 

y[t] = a(Jy[t-l] + g)- (39) 

with the starting condition y[0] = 0. Results are shown for numbers of iterations ranging from 1 to 
19. As few as 5 iterations gives a good approximation for small c, while increasing the number of 
iterations to 19 adds little. The low-contrast iterations all fail before or very slightly after c = 10, 
which corresponds to a in the range 1.4 to 2.4 across the parameters. That is, the failure occurs 
for a = 0(1), as expected. 

For the high-a or small-/? (high contrast) case, we treat the equation y = ^— g -|- fiy'n^ as 
the recurrence relation: 

y[t]=J-i(-g + ;3y[t-l]-^) (40) 

We use as starting conditions ?/e[0] = with yi{0\ = for Q^; > 0, y/[0] = qi/Jei for < 0. 
For rig < 0, using y[0] = would give complex solutions. We instead use as a starting condition 
the value of y when yE reaches zero with increasing c. The fourth row of Fig. 2 illustrates these 
high-contrast solutions. Again, 5 iterations do about as well as larger numbers of iterations. The 
iterations give good approximations for high c but, for VLe < 0, fail for larger c as te approaches 
zero. (3 is very small for these large c's, so this presumably represents the initial conditions no longer 
being in the basin of attraction of the fixed point. For low c failure of convergence is expected for 
P = 0(1) (recall that /3 increases with decreasing c, /3 = 1/y/a), although problems with the basin 
of attraction could also arise. None of the iterations work for c below about 9 or 10, corresponding 
to /3 roughly in the range .65 to .85, with the exception of column E. In that column, the largest 
number of iterations works down to the jump in ve and rj, which occurs at about c = 5.435 for 
the given parameters, or f3 around 1.1. In column D the iterations do not work below c about 190, 
which corresponds to (3 about 0.15, a bit lower than expected. 

Finally, in the fifth row of Fig. 2, we show the value of q = ^ that divides stability (values 
below curves) from instability (values above curves) of the fixed point, according to Eq. 29. In all 
cases except $7/ < il^; < 0, the fixed point remains stable for q < 1 across the range of studied 
stimulus strengths, indicating that fine tuning or unreasonably small values of q are not required. 

6 Discussion 

We have shown in studies of a 2-D system (and found in simulation studies of higher-dimensional 
systems, to be presented elsewhere) that the supralinear network will dynamically stabilize with 
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increasing input strength provided the I ^ E and E ^ I connections mediating feedback inhibition 
arc sufficiently strong and the inhibitory time constant is not too slow. This dynamic stabilization 
results in a change from responses scaling supralinearly to responses scaling sublinearly with input 
strength. The system can also yield "supersaturation" , in which excitatory firing rates reach a peak 
with increasing input strengths and then decrease (as observed biologically, Ledgeway et al. 2005, 
Li and Creutzfeldt 1984, Peirce 2007, Tyler and Apkarian 1985), with rates ultimately decreasing 
to zero for large enough input strengths (which presumably are beyond the dynamical range of 
biological inputs). The conditions for this to occur were characterized in the 2-D system. The 
strongest sublinear behavior, and hence behavior most likely to underly biological observations in 
cerebral cortex, occurs for parameters that lead to supersaturation. 

Many questions remain outstanding. As some examples: within the range of models analyzed 
here, can more precise results, analogous to those obtained here for 2-dimensional models, be 
obtained for higher-dimensional models, for which we only discussed general scaling arguments? 
For any dimensionality, can useful results be obtained as to when the network is globally stable? 
How will diversity of network parameters, including in particular of the power n, alter behavior? 
Presumably an even slightly larger mean n for I vs. E cells will enormously enhance the range of 
parameters that will stabilize. How will cell-to-cell variability of n affect behavior? How will behav- 
ior be affected by taking into account the decreased noise level, and thus increase in n, that occurs 
with increasing stimulus contrast (Finn et al. 2007), i.e. with increasing input firing rate? How will 
network behavior be modified by addition of short-term synaptic facilitation and depression {e.g., 
Fioravante and Regehr 2011)? Can analysis be done of more biophysically realistic models, such 
as networks of integrate-and-fire neurons, which have an input/output function well approximated 
by a power law so long as they are firing on input fiuctuations rather than the mean input (Hansel 
and van Vreeswijk 2002)? What can we learn as we move beyond the steady state to network 
dynamics, particularly using more realistic models that can better capture faster dynamics and 
that incorporate synaptic delays? How will the network behave when multiple types of inhibitory 
neurons {e.g. Isaacson and Scanziani 2011), or of excitatory neurons, are incorporated? Incorpo- 
rating multiple neuronal subtypes must presumably be guided by knowledge not yet available of 
the separate connectivity patterns as well as biophysical properties of the different subtypes. 

Despite the many open questions, we believe the basic findings are likely to be quite robust 
and to underly a wide range of cerebral cortical behavior: networks with supralinear input/output 
functions can dynamically stabilize, resulting in a transition from supralinear to sublinear input 
summation. 
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